function [pv,pIv,vv,jstar,bfail] =SuffNSolBellf(PI,M,vbarv,bet,gam) 
% Bellman function iteration N-state model

N = max(size(vbarv));
bfail=1;
pbarv=(eye(N)-bet*gam*PI)\vbarv;


pI = ((PI.*M)*gam)*pbarv;
pbarU = pbarv;
palg = max(pI,pbarU);

k=200000;

critm1 = 10e-9;
%Palg=zeros(N,k);
for j=1:k
    pI = ((PI.*M)*gam)*palg;
    pbarU =vbarv+ bet*gam*PI*palg;    
    palgn = max(pI,pbarU);
    if j>500
        diffm1=max(abs(palgn-palg));
        if diffm1<critm1
            'Bellman equation successfully converged'
            j
            palg=palgn;
            bfail=0;
            break
        end
    end
    palg=palgn;
end

pv=palg;
pIv = pI;
vv = max(vbarv,(((PI.*M)*gam)-(bet*gam)*PI)*pv);
jstar = sum(pv-pIv==0);  %last state in -- 1 factor model

end

